# -*- coding: utf-8 -*-
from datetime import datetime, time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.metrics import mean_squared_error
from numpy import sin, exp, cos

plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False

def plot_month_data(data):
    fig, ax = plt.subplots()
    ax.plot(data)
    ax.set_xlabel('日期', fontsize=20)
    ax.set_ylabel('接警次数', fontsize=20)
    
    ax.legend()
    ax.grid()

    fig.savefig(r"../图片/1.png")
    plt.show()
    
def plot_evaluation(data, data_test, model_name):
    '''
    画出拟合数据和测试数据
    '''
    fig = plt.figure()
    plt.plot(np.arange(49,61), data, label='Fitted Data')
    plt.plot(np.arange(49,61), data_test, label='Real Data')
    plt.xlabel('日期序列', fontsize=16)
    plt.ylabel(model_name, fontsize=16)
    plt.grid()
    plt.legend()
    fig.savefig(f"../图片/{model_name}.png")
    plt.show()
    
def linear(x):
    '''
    一元线性回归
    '''
    y = 0.01485*x + 62.16
    return np.round(y)

def polynomial(x):
    '''
    多元线性回归
    '''
    y = -0.006076*x**3 + 0.418*x**2 - 7.473*x + 88.15
    return np.round(y)

def sum_sin(x):
    '''
    sin 函数累加
    '''
    a1 =       262.4  
    b1 =    0.005264  
    c1 =       2.781  
    a2 =       28.32  
    b2 =      0.1117  
    c2 =        3.83  
    a3 =       22.68  
    b3 =      0.5016  
    c3 =      0.3463  
    a4 =       10.87  
    b4 =      0.3486  
    c4 =     -0.1735 
    a5 =       12.85 
    b5 =       1.586 
    c5 =      0.2025 
    a6 =       16.22  
    b6 =       1.043
    c6 =       2.913  
    a7 =        16.5  
    b7 =       2.072  
    c7 =     -0.9097  
    a8 =       9.579  
    b8 =       1.483  
    c8 =      -3.065  

    y = a1*sin(b1*x+c1) + a2*sin(b2*x+c2) + a3*sin(b3*x+c3) + \
        a4*sin(b4*x+c4) + a5*sin(b5*x+c5) + a6*sin(b6*x+c6) + \
        a7*sin(b7*x+c7) + a8*sin(b8*x+c8)

    return np.round(y) 
    
def fourier(x):
    '''
    傅里叶函数
    '''
    
    a0 =       63.39  #(57.41, 69.37)
    a1 =       6.254  #(-2.21, 14.72)
    b1 =       2.553  #(-5.914, 11.02)
    a2 =     0.03952  #(-9.317, 9.396)
    b2 =          26  #(17.5, 34.5)
    a3 =       8.743  #(0.1131, 17.37)
    b3 =       4.235  #(-4.65, 13.12)
    a4 =      0.6468  #(-9.578, 10.87)
    b4 =      -13.69  #(-22.2, -5.166)
    a5 =      -3.247  #(-12.57, 6.078)
    b5 =       9.375  #(0.449, 18.3)
    a6 =       16.76  #(8.055, 25.46)
    b6 =       5.683  #(-6.345, 17.71)
    a7 =      -4.149  #(-13.34, 5.04)
    b7 =      -1.762  #(-11.05, 7.526)
    a8 =      -10.56  #(-21.23, 0.1067)
    b8 =       10.92  #(-0.6063, 22.45)
    w =      0.2582  #(0.2547, 0.2618)

    y = a0 + a1*cos(x*w) + b1*sin(x*w) + \
       a2*cos(2*x*w) + b2*sin(2*x*w) + a3*cos(3*x*w) + b3*sin(3*x*w) + \
       a4*cos(4*x*w) + b4*sin(4*x*w) + a5*cos(5*x*w) + b5*sin(5*x*w) + \
       a6*cos(6*x*w) + b6*sin(6*x*w) + a7*cos(7*x*w) + b7*sin(7*x*w) + \
       a8*cos(8*x*w) + b8*sin(8*x*w)
       
    return np.round(y)

def gaussian(x):
    '''
    高斯函数拟合
    '''
    a1 =         372  #(-1.403e+006, 1.404e+006)
    b1 =       16.44  #(-136.1, 169)
    c1 =      0.4019  #(-497, 497.8)
    a2 =       127.8  #(96.3, 159.3)
    b2 =       40.64  #(40.43, 40.86)
    c2 =       1.106  #(0.7444, 1.468)
    a3 =       108.7  #(-2323, 2540)
    b3 =       25.29  #(18.19, 32.39)
    c3 =      0.5491  #(-8.553, 9.652)
    a4 =       76.62  #(46.67, 106.6)
    b4 =       28.51  #(28.14, 28.88)
    c4 =       1.261  #(0.6203, 1.903)
    a5 =       74.79  #(47.94, 101.6)
    b5 =      0.5594  #(-3.994, 5.113)
    c5 =       6.917  #(1.192, 12.64)
    a6 =       35.59  #(9.232, 61.95)
    b6 =       36.85  #(35.9, 37.81)
    c6 =       1.651  #(0.07488, 3.227)
    a7 =       27.15  #(-9.402, 63.71)
    b7 =       22.74  #(20.35, 25.12)
    c7 =        1.17  #(-3.037, 5.377)
    a8 =       47.41  #(31.71, 63.12)
    b8 =       29.44  #(23.6, 35.27)
    c8 =       22.96  #(9.203, 36.71)
    
    y = a1*exp(-((x-b1)/c1)**2) + a2*exp(-((x-b2)/c2)**2) + \
              a3*exp(-((x-b3)/c3)**2) + a4*exp(-((x-b4)/c4)**2) + \
              a5*exp(-((x-b5)/c5)**2) + a6*exp(-((x-b6)/c6)**2) + \
              a7*exp(-((x-b7)/c7)**2) + a8*exp(-((x-b8)/c8)**2)
              
    return np.round(y)
    
def evaluate_model(data_test):
    x = np.arange(49, 61, 1)
    
    y_lr = linear(x)
    y_poly = polynomial(x)
    y_sin = sum_sin(x)
    y_fourier = fourier(x)
    y_gaussian = gaussian(x)
    
    model = {'y_lr':'线性模型', 'y_poly':'多元模型', 'y_sin':'sin函数累加模型',\
            'y_fourier':'傅里叶模型', 'y_gaussian':'高斯函数模型'}
    
    for var_name, model_name in model.items():
        
        data_fit = vars()[var_name]
        score = mean_squared_error(data_fit, data_test)
        plot_evaluation(data_fit, data_test, model_name)
        print(f'{model_name} 的均方误差为： {score}')
    
    
if __name__ == '__main__':
    path = r'../附件/附件2：某地消防救援出警数据.xlsx'
    time_series = pd.read_excel(path)
    # 根据月份找出事故次数
    data = time_series.groupby(pd.Grouper(freq='1M', \
                                key='接警日期')).count().iloc[:, 1].values
    
    plot_month_data(data)
    
    # 将数据导出为 excel
    month_data = pd.DataFrame(data=data, columns=['事故次数'])
    month_data.to_excel(r'../附件/事故次数(月).xlsx')
    
    # 判断各一元模型的效果
    month_data_test = month_data[48:].values.reshape(12,)
    
    # 评价各个一元模型的效果，评价方法为 MSE：
    
    evaluate_model(month_data_test)
    
